% ML TAU 2013 final project script
% Grid search:
% Calculate abs(S1 - S2) and perform SVM based on its L to H PCA

function grid_search2()

%
% Initialization
%
close all; clear; clc; tic;

%base_path = '/home/itay/TAU/IML/final';
libsvm_path = './libsvm-3.17/matlab';

% Add path to the libsvm
%addpath(base_path);
addpath(libsvm_path);

N = 21; %number of results to be saved.

if ~exist('dataforproject.mat','file')
    error('Data file not found in current directory')
end;

% Saving memory. Loading vars on a need to load basis only throughout
load('dataforproject.mat','X1train','X2train','gidtrain','ytrain');
fprintf('Training Data successfully loaded\n');

X_DELTA = abs(X1train-X2train);

L_vec = 3:3; %2:5

H_vec = 120:120; % [50,80,100,120,150,200,500,1000];

ker_degrees = 1:6;

ker_types = [3,1];

c_degrees = [-4,14];

best_results = zeros(N,6) ; i = 1;

[L, H, C, D, T, Err] = svm_grid_search(ker_types,ker_degrees,c_degrees,L_vec,H_vec);
fprintf('===\tResults\t===\n(L*=%d ; H*=%d ; C*=%f ; t*=%d ; deg*=%d)\n#Avg Acc: %f\n',L,H,C,T,D,1-Err);
save(sprintf('best_%d_results.mat',N), 'best_results');
disp(best_results);

    function [L, H, C, D, T, Err] = svm_grid_search(ker_types,ker_degrees,c_degrees,L_vec,H_vec)
        
        s = size (c_degrees);
        
        c_vec = power(repmat(2,s),c_degrees);
        
        Err = 1.0; L = L_vec(1); H = H_vec(1);
        
        for l=L_vec
            
            for h=H_vec
                
                
                % Calculate and cache the test and train data for cross validation
                
                num_train_entries = size(X_DELTA(:, (gidtrain~=1)), 2);
                num_test_entries = size(X_DELTA(:, (gidtrain==1)), 2);
                
                num_pcs = (h - l) + 1;
                NormalizedTrainX = zeros(num_train_entries, num_pcs, 3);
                NormalizedTestX = zeros(num_test_entries, num_pcs, 3);
                y_train = zeros(3, num_train_entries, 1)';
                y_test = zeros(3, num_test_entries, 1)';
                
                for k=1:3
                    % Training set
                    X = X_DELTA(:, (gidtrain~=k));
                    y_train(:, k) = ytrain(gidtrain~=k);
                    
                    [ EigenVectors, MeanX ] = PerformPCA(X, [ 'X_DELTA', int2str(k) ]);
                    TransformedX = ((X - repmat(MeanX, 1, size(X, 2)))' * EigenVectors);
                    TransformedX = TransformedX(:, L:H);
                    
                    % Scale X to [ -1, 1 ]
                    Minimums = min(TransformedX, [], 1);
                    Ranges = max(TransformedX, [], 1) - Minimums;
                    NormalizedTrainX(:, :, k) = (2 * (TransformedX - repmat(Minimums, num_train_entries, 1)) ./ repmat(Ranges, num_train_entries, 1)) - 1;
                    
                    % Testing set
                    X = X_DELTA(:, (gidtrain==k));
                    TransformedX = ((X - repmat(MeanX, 1, size(X, 2)))' * EigenVectors);
                    TransformedX = TransformedX(:, L:H);
                    y_test(:, k) = ytrain(gidtrain==k);
                    
                    % Scale to [ -1, 1 ]
                    NormalizedTestX(:, :, k) = (2 * (TransformedX - repmat(Minimums, num_test_entries, 1)) ./ repmat(Ranges, num_test_entries, 1)) - 1;
                end
                
                for  t=ker_types
                    
                    for d=ker_degrees
                        
                        for c=c_vec
                            
                            Results = zeros(3, 1);
                            
                            for k=1:3
                                
                                model = svmtrain(y_train(:, k),NormalizedTrainX(:, :, k), sprintf('-t %d -g 1.0 -c %f -d %d -h 0 -q', t, c, d));
                                
                                % Predict
                                predictedY = svmpredict(y_test(:, k),NormalizedTestX(:, :, k),model);
                                
                                Results(k) = sum(y_test(:, k).*predictedY > 0) / size(y_test(:, k),1);
                                
                            end
                            
                            accuracy = mean(Results);
                            fprintf('Accuracy for C=%f D=%d T=%d H=%d L=%d: %f\n', c, d, t, h, l, accuracy);
                            
                            e = 1-accuracy;
                            if e < Err
                                C=c; D=d; T=t; Err=e; H=h; L=l;
                                fprintf('**** New best error rate found: C=%f D=%d T=%d Err=%f H=%d L=%d\n', c, d, t, e, h, l);
                            end
                            if e < 0.3
                                best_results(mod(i,N)+1,:) =[l h c d t e];
                                i = i+1;
                            end
                        end
                    end
                end
            end
        end
        
    end

end

